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ABSTRACT 

A technique for the construction of axisymmetric distribution functions for individual 
galaxies is presented. It starts from the observed surface brightness distribution, which 
is deprojected to gain the axisymmetric luminosity density, from which follows the 
stars' gravitational potential. After adding dark mass components, such as a central 
black hole, the two-integral distribution function (21- DF) f(E, L z ), which depends only 
on the classical integrals of motion in an axisymmetric potential, is constructed using 
the Richardson-Lucy algorithm. This algorithm proved to be very efficient in finding 
f(E, L z ) provided the integral equation to be solved has been properly modified. Once 
the 2I-DF is constructed, its kinematics can be computed and compared with those 
observed. Many discrepancies may be remedied by altering the assumed inclination 
angle, mass-to-light ratio, dark components, and odd part of the 2I-DF. Remaining 
discrepancies may indicate, that the distribution function depends on the non-classical 
third integral, or is non-axisymmetric. 

The method has been applied to the nearby elliptical galaxy M32. A 2FDF with 
~ 55° inclination and a central black hole (or other compact dark mass inside ~lpc) 
of 1.6-2 x 1O 6 M fits the high-spatial-resolution kinematic data of van der Marel et al. 
remarkably well. 2I-DFs with a significantly less or more massive central dark mass or 
with edge-on inclination can be ruled out for M32. Predictions are made for observations 
with the HST: spectroscopy using its smallest square aperture of 0'.'09 x 0'.'09 should 
yield a non-gaussian central velocity profile with broad wings, true and gaussian-fit 
velocity dispersion of 150-170 kms -1 and 120-130 kms -1 , respectively. 

Key words: stellar dynamics - galaxies: kinematics and dynamics - galaxies: structure 
- galaxies: central black holes - galaxies: individual (M32) - line: profiles - methods: 
numerical 
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1 INTRODUCTION 



Since the 1970s it is known that, in general, elliptical gala- 
xies are not simply isotropic oblate stellar systems flattened 
by rotation, but have triaxial shapes, which are supported 
by anisotropic velocity distributions. To better understand 
the dynamical structure of these objects, detailed dynami- 
cal models are required, and a natural way of construct- 
ing them consists of the following steps, (i) Find a three- 
dimensional luminosity density p which gives the observed 
surface brightness when projected and convolved with the 
atmospheric seeing; (ii) assume a mass-to-light ratio and 
integrate the Poisson equation to obtain the potential of 
the stars, to which further mass components may be added, 
e.g. a central black hole or a dark halo; (iii) solve the in- 
tegral equation p = J fd 3 v for the distribution function /, 
which must depend only on the isolating integrals of motion 
in the given potential to satisfy the collisionless Boltzmann 
equation; and (iv) finally, compare the kinematic observables 



with those observed, to test whether the model is consistent 
or not. 

In general, this procedure is very complicated and al- 
lows a huge amount of freedom. Already the deprojection of 
a triaxial body is highly non-unique. However, if one assumes 
axisymmetry and fixes the inclination, there is de facto only 
one smooth and physical density distribution for any allowed 
inclination (Palmer 1994), but strictly speaking the depro- 
jection is unique only for edge-on inclination (Rybicki 1987). 
Another reason for the assumption of an axisymmetric con- 
figuration is its simple phase-space structure easing the dy- 
namical modelling: there is only one major orbit family and 
two classical integrals of motion, while triaxial potentials al- 
low for four major orbit families with, in general, the energy 
being the only classical invariant. Additionally, there is also 
some observational evidence, that elliptical galaxies in their 
majority have near-oblate shapes, as indicated by the align- 
ment of kinematic and photometric axes in many of these 
objects (Franx, Illingworth & de Zeeuw 1991). Another in- 
dication comes from studies of gas flows, for instance, in the 
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case of the elliptical galaxy IC 2006 the form and velocity 
field of an equatorial gas ring constrains the potential to 
be round in the equatorial plane (Franx, van Gorkom & de 
Zeeuw 1994). This view is supported by recent results of nu- 
merical iV-body simulations with small amounts of gas being 
included: they predict near-oblate shapes both for remnants 
of disk mergers (Barnes & Hernquist 1994) and for dark mat- 
ter halos in CDM cosmogonies (e.g. Katz & Gunn 1991). 

Typical stellar orbits in axisymmetric potentials obey 
three integrals of motion, namely the energy E, the z- 
component of angular momentum L z , and the non-classical, 
third integral /3. Unfortunately, the third integral is gen- 
erally not known explicitly and must be evaluated by per- 
turbation theory or similar means. Therefore, one often re- 
stricts oneself to models with stellar distribution functions 
(DFs) depending only on the two classical integrals of mo- 
tion, i.e. / = f(E,L z ) (2I-DF). Though there is no obvi- 
ous physical reason for the DFs of axisymmetric galaxies 
to ignore the third integral, such models are useful for the 
following reasons, (i) The part of the DF even in L z , f e , fol- 
lows uniquely from density and potential, and hence from 
the observed surface brightness and the assumptions made 
on inclination, mass-to-light ratio, and dark mass compo- 
nents. This substantially reduces the amount of freedom in 
the DF. (ii) Since a two-integral model (2I-model) is simply 
to construct (compared to any 3I-model), it can be used as a 
template, i.e. the comparison of a galaxy's kinematics with 
that predicted by the 2I-model already gives clear hints as 
to the true nature of the dynamical configuration of that ga- 
laxy. Instead of modelling the DF, one may also consider its 
velocity moments, i.e. streaming velocity and velocity dis- 
persion, which in the case of / = f(E, L z ) follow uniquely 
from the Jeans equations. However, there is always the pos- 
sibility that the underlying DF is not positive definite, which 
may be undetectable from its moments. Moreover, nowadays 
not only velocity and dispersion of the stellar line-of-sight 
velocity profiles (VPs), but also their shapes can be extracted 
from galaxy spectra (Bender 1990; Rix & White 1992; van 
der Marel & Franx 1993; Winsall & Freeman 1993; Kuijken 
& Merrifield 1993). To fully exploit these data one must ei- 
ther calculate the DF or solve the Jeans equations for the 
higher order (~16) velocity moments (Magorrian & Binney 
1994). 

In this paper I present a numerical method for the con- 
struction of two-integral distribution functions for individual 
galaxies. These allow to explore the significance of specific 
observational features for the dynamical structure of a ga- 
laxy. Thus such models grant deeper insights than do sim- 
ple spheroidal models as those of Dehnen & Gerhard (1994) 
or Qian et al. (1994). The method essentially consists of 
(1) a deprojection of the surface brightness into an axisym- 
metric luminosity density and (2) the construction of the 
corresponding 2I-DF with the assumption of some mass-to- 
light ratio, which may include dark mass components such 
as a massive central black hole or a dark halo. Both steps 
are similar insofar as a linear integral equation has to be 
solved. In the given method this is done using the Richard- 
son (1972) -Lucy (1974) (RL) algorithm, which - after some 
modifications - turns out to be very efficient in solving both 
equations. 

The method is illustrated by an application to the 
nearby elliptical M32, which is suspected harbouring a mas- 



sive central black hole (Tonry 1987; Dressier & Richstone 
1988; Richstone, Bower & Dressier 1990). Recently, Qian 
et al. (1994) have employed the method of Hunter & Qian 
(1993) to compute f(E, L z ) for a simple mass model for the 
central ~20" of M32 used by van der Marel et al. (1994b) in 
modelling the velocity moments. Although the model pre- 
sented here differs from the Qian et al. analysis in details 
of the stellar mass distribution and in the technique for the 
evaluation of f(E,L z ), the main result is the same: a 21- 
model with a 1.6-2 x 10 6 Mq black hole gives a good fit to 
the kinematic data observed by van der Marel et al. (1994a). 
However, some small discrepancies between the VPs of M32 
and of the models remain, which suggest that the DF may 
depend weakly on the third integral. 

The sections of this paper are designed to be read al- 
most independently from each other. Their contents are as 
follows. Section 2 describes how the RL algorithm can be 
used to find f(E, L z ) from density and potential; numerical 
details and the result of a test with an analytical model are 
given. The construction of two-integral models for individual 
galaxies and the comparison of their kinematics with that 
observed is outlined in Section 3. As an application Section 
4 deals with 2I-models for M32, and Section 5 sums up and 
concludes. Appendix A briefly describes the RL algorithm 
and how to achieve faster convergence; Appendix B contains 
formulae for the recovery of f(E, L z ) from the derivatives of 
the density; and Appendix C gives a Monte-Carlo procedure 
for the evaluation of VPs accounting for binning and seeing. 



2 NUMERICALLY RECOVERING f(E, L z ) 
FROM DENSITY AND POTENTIAL 

For an axisymmetric stellar system with a 2I-DF the spatial 
density p is related to the phase space density / by 

P (*,R) = 2nJ AE J ^ fe{E ,L 2 z ), (1) 

where the density is written as a function of the cylindri- 
cal radius R and the relative potential , F, while f e denotes 
the part of the 21- DF even in L z , which only contributes to 
the density. It has been known for more than 30 years that 
the solution of this equation for f e (E,L 2 z ) given p(JB, R) is 
unique (Lynden-Bell 1962), however, only recently has a ge- 
neral inversion formula been found (Hunter & Qian 1993), 
that is the generalization of Eddingtons (1916) formula for 
the f(E) generating a spherical system. This inversion for- 
mula requires the knowledge of dp/d^fr in parts of the com- 
plex 'F-plane, which can be numerically evaluated for many 
functional forms of p, even if it is not expressed by elemen- 
tary functions in 'F and R (Qian et al. 1994). However, this 
technique can hardly be applied to densities in tabulated 
form, since (i) the evaluation of dp/d^fr is almost impossible 
with high accuracy, and (ii) the density may well be contam- 
inated with noise, which is amplified to strong oscillations 
in the DF. 

In this section a scheme to solve for f e (E,L 2 z ) will be 
given that is designed for tabulated density fields, guarantees 
a physical DF, and nevertheless is fast and accurate. 
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2.1 Using the Richardson-Lucy algorithm 

A suitable tool for the recovery of f e (E, L 2 ) for general 
axisymmetric density fields is the Richardson (1972) -Lucy 
(1974) algorithm (hereafter RL algorithm). It was first in- 
troduced into stellar dynamics by Newton & Binney (1987), 
who used it to construct a positive spherical DF for M87; 
later on Gerhard (1991) applied it for the evaluation of an- 
isotropic DFs in spherical potentials. A description of this al- 
gorithm, and how to accelerate and modify it for faster con- 
vergence, is given in Appendix A. Starting from any guess 
for f e (E, L 2 ), one first projects the DF into p(JS , R) accord- 
ing to equation (1), and then corrects f e (E, L 2 ) using the 
discrepancy between the actual density and that of the DF. 
This pair of steps may then be repeated until satisfactory 
convergence is obtained. Applying the RL algorithm after 
multiplying equation (1) with R~ s p~ p on both sides (Ap- 
pendix A2) yields for the correction step 



fe,n+l (E, L 2 Z ) 
fe,n(E, L 2 Z ) 



* 2 rR 
d* 
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The numbers s and p may be chosen to optimize conver- 
gence; Ry, is given by ^l/(i?$,0) = , F; $i and $2 are the 
roots of * = E + L 2 Z /(2R%); and p n (^,R) is the density 
generated by f etn (E, L 2 ), the DF after the ra-th iteration 
step. For L z = L c i vc (E) the double integral shrinks to that 
point to which the circular orbit contributes. For L z = the 
integral over R diverges at R = 0. However, one may apply 
the RL algorithm to p(#, 0) = 4tt /* d£^2(* - E)f(E, 0) 
to get 
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Thus, recovering f e (E,L 2 z ) from p(JS , R) via the RL algo- 
rithm yields a procedure, which at each iterative step re- 
quires a two-dimensional integration on almost each point of 
the two-dimensional grids in (E,L 2 Z ) and (JS,R). Appendix 
B describes procedures, which require only one-dimensional 
integrations at each point of the two-dimensional grids. How- 
ever these involve taking the derivative with respect to 'F or 
R on both sides of equation (1), while my goal is to recover 
f e (E,L 2 z ) from a tabulated density, for which these deriva- 
tives cannot be evaluated with sufficient accuracy. 

2.2 Unphysical distribution functions 

There may well be combinations of positive density and po- 
tential for which no underlying non-negative f e (E,L 2 z ) ex- 
ists. Then the scheme given above will never converge to 
p = p n (within the numerical accuracy) as n increases. The 
remaining discrepancies in the density may either by con- 
fined to spatial scales as small as the resolution element of 
the grid used to represent the density, or they occur on a 
larger scale. In the first case, the discrepancies may easily 
be interpreted as noise, since the RL algorithm guarantees 
a positive DF, and hence gives a smooth density. 

In the second case, the large scale behaviour of density 
and potential obliges the true f e (E,L 2 z ) to be somewhere 



negative. When using the RL scheme, this will be recognized 
by the fact, that the error in the density does not converge 
to zero, despite the DF continuously decreases in a certain 
region of phase space. In such a case, the above scheme can 
be used to recover the true but non-physical f e (E,L 2 z ) by 
applying it a second time on the density p = kp n — p, where 
k has to be chosen to keep p everywhere positive, and fi- 
nally constructing the non-physical DF as difference of two 
physical ones. 



2.3 Numerical implementation 

For the numerical implementation of the procedure, p(^,R) 
is tabulated on a grid logarithmic in Ry, and equidistant in 
R/R<$, whereas f e (E, L 2 ) is tabulated on a grid logarithmic 
in Re, where ^(Re, 0) = E, and equidistant in L 2 /L 2 ilc (E). 
Interpolations are done in Inp and In/ using cubic splines. 
The integrals dL 2 and dR in equations (1) and (2) are com- 
puted by integrating the spline curves exactly, while those 
over E and 'F are evaluated by gaussian quadratures. 

Starting from a density field p(R, z) and a potential 
^(R, z), the routine first eliminates z to yield p(^, R), com- 
putes $1 and ^2 for the above grid in (E,L 2 Z ), and pre- 
calculates the integrals in the denominators of equations (2) 
and (3) on that grid. An initial f(E, 0) is guessed using 
Eddingtons (1916) formula for isotropic spherical systems, 
which is also valid for f(E, 0) of axisymmetric systems. The 
initial guess for f e (E, L 2 ) is then obtained by first using the 
RL algorithm to fully recover f(E, 0) (equation 3) and then 
multiplying it with some function of L 2 / L 2 1TC . The conver- 
gence of the iteration is accelerated as described in Appendix 
A3. After each correction step, / e , n is smoothed on the scale 
of the grid cells, to avoid fitting small-scale noise. With a 
grid of 80x11 points for the density, 131x17 points for the 
DF, and 60 points for the gaussian quadratures, each iter- 
ation takes about 40 seconds CPU time on a DEC Alpha 
3000-600 workstation. 



2.4 Testing with a simple model 

Recently, Evans (1994) found a family of axisymmetric mo- 
dels with simple analytical 2I-DFs. In these models the po- 
tential is proportional to (R 2 + R 2 + z 2 q~ 2 )~^^ 2 , and the 
21- DFs are always finite sums of powers in E and L 2 Z . To test 
the algorithm and its numerical implementation, I chose the 
model with fi = 0.5 and q = 0.93. The numerical apparatus 
was fed by the density field tabulated between 0.005i? c and 
200i? c , the potential was obtained by multipole expansion. 
After only three iterations the routine found f e (E, L 2 ), such 
that the RMS relative deviation between its density and the 
input was as small as 0.003. Figure 1 compares the exact 
DF (solid contours) with the numerical results (dashed), the 
differences are almost invisible. 



3 MODELLING INDIVIDUAL GALAXIES 
WITH TWO-INTEGRAL MODELS 

This section describes how two-integral models for individ- 
ual galaxies can be obtained and compared to the observa- 
tions. A first application follows in Section 4 with the E3 
galaxy M32. 
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Figure 1. Contour lines of f e (E,L 2 ) for Evans' (1994) power-law 
model with f3 = 0.5 and q = 0.93; solid: analytical DF, dotted: DF 
numerically reconstructed via the RL algorithm. 



I(x, y) 



4groPo 



I 2 . 2 — 2 \ / 2 

(x +y q P )/r , 



1 ^ + t 2 y-'(i-t 2 f- 2 dt 



2 

<lp 



^ 2 f + (1 - f )* 2 

2 • 2 • . 2 • 

: g sin i + cos i. 



(7) 



The parameters 7, fi , po and ro are fitted to the surface 
density along the major axis, the axis ratio q is guessed 
from a mean value for the apparent axis ratio q p and the 
inclination. 

In order to ensure, that the errors are not systemati- 
cally distributed, i.e. that the result is not biased towards 
the starting density, a minimum number of iterations is nec- 
essary. To avoid an amplification of noise in the data, the 
density field may be smoothed after each correction step on 
a scale similar to the resolution element of the grid. 

Finally the stellar mass distribution is obtained by mul- 
tiplying the luminosity density with a (usually constant) 
mass-to-light ratio T, and the stars gravitational potential 
can be evaluated, e.g. by using multipole expansion (cf. Bin- 
ney & Tremaine 1987, §2.4). Additional dark mass compo- 
nents may be added to the potential. 



3.1 The mass model 

There are two fundamentally distinct ways to derive a three- 
dimensional luminosity density from a given surface bright- 
ness: fitting some parameterized function and using unpa- 
rameterized methods. Recently, Merritt & Tremblay (1994) 
have discussed the relative merits of both techniques. In the 
present framework, a fitted parameterized function would be 
free of small scale noise and advantageous in computing the 
gravitational potential, however, the density would always 
be biased. To allow for general forms of p(R, z) (i.e. any ra- 
dial run, disky/boxy contours, changing ellipticity), I use an 
(almost) unbiased technique: the RL algorithm, which was 
introduced for this purpose by Binney, Davies & Illingworth 
(1990). Differently from these authors, I modify the correc- 
tion equation with a = I~ p (in the notation of Appendix 
A2), to yield a better convergence. Then the projection and 
correction equation are 

I„(x,y) = — — : / ds p n (\/ s 2 + x 2 , y seci + s coti) (4) 
smi / 

•J — CO 

in _\ f 7 ^ 2 , dt I„(Rcost, z smi + R cosi smt) 

Pn-\-l l^Jl-, ZJ j_7t/2 v ' 



Pn (R,z) r«/ 2 



J_ n , 2 dt[I(R cost, z smi + R cosi smt)] p 



(5) 



with I n (x,y) = I(x, y) 1 ~ p I I n (x, y), where x and y are re- 
spectively the co-ordinates in the plane of sky along the ma- 
jor and minor axes of the galaxy. The angle of inclination is 
denoted i with 90° meaning edge-on projection. The choice 
p = 0, for which the above equations reduce to those already 
given by Binney et al., turned out to give best convergence 
in the case of a centrally flat surface brightness, while for a 
cusp p fa 1 is more efficient. As initial guess for the density 
I use the spheroidal model: 



3.2 The velocity profiles 

After the construction of f e , as outlined in Section 2, one 
finally can compute its observables. The complete set of 
kinematic observables for a stellar system are the velocity 
profiles (hereafter VP) at each position on the sky. The VP 
l(vi os ;x,y) is simply the phase space distribution function 
integrated over the non-observable co-ordinates, which are 
the stars' positions s along the line of sight and their veloc- 
ities in the plane of the sky 



l(v los ;x, y) 



ds I dv x I dv y f(E, L z ). 



(8) 



p(R, z) = p m 7 (m + l) 



7-/3. 



(R 2 + z 2 q~ 



! )Ao, (6) 



Robust algorithms to extract the VPs from observed galaxy 
spectra have been available for a few years (Bender 1990; 
Rix & White 1992; van der Marel & Franx 1993; Kuijken 
& Merrifield 1993; Winsall & Freeman 1993). In general, 
spectra are only measured at a few slit positions, i.e. major 
and minor axes, and only in the central regions of galaxies 
are the signal-to-noise ratios high enough to extract the VPs. 
In order to compare the model with such observations, one 
must account for effects of the finite slit width and pixel 
size as well as for atmospheric seeing. All these will lead to 
some averaging and can change the VPs near the center of a 
galaxy dramatically. Appendix C gives a numerical method 
for the computation of VPs of 2I-models based on Monte 
Carlo methods, which allow these effects to be included. 



3.3 Comparing with the observed kinematics 

Because of the uniqueness of f e , it is useful, to first compare 
its observables with those of the modelled galaxy and if nec- 
essary correct the inclination or mass-to-light ratio before 
calculating f , the DF's odd part. The observables of f e are 
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simply the even parts of the VPs, VP e , and its moments 1 
especially the RMS-velocity {v 2 ^) 1 ^ 2 . 

The RMS-velocity scales with the square root of the 
(constant) mass-to-light ratio T, while the inclination in- 
fluences the ratio of RMS-velocities on the minor and major 
axes ("i 2 os )min/("i 2 os}maj- The la tter decreases with decreas- 
ing i<90° (Binney et al. 1990; Dehnen & Gerhard 1994). 
A central black hole alters the kinematics in the central re- 
gions by causing {v 2 ^) 1 ^ 2 to rise as r -1 ^ 2 (neglecting seeing), 
while a dark halo alters the radial run of {v 2 ^) 1 ^ 2 at large 
radii. The behaviour with inclination implies, that a galaxy 
with a greater ratio ("i 2 os )min/("i 2 os}maj than a 2I-model with 
assumed edge-on projection cannot have a 21- DF, and there- 
fore must be described by a three-integral model 2 . 

After inclination, mass-to-light ratio and any dark mass 
components have been adapted to give {v 2 os Y^ 2 in agreement 
with the observations, the shapes of the VP e can be com- 
pared with the observed ones, which again may rule out a 
21- DF. Once f e has been determined, the remaining job is to 
test, whether there exists a f that (i) leads to a physical DF 
and (ii) results in VPs in agreement with the observations. 
One constraint is that the observed mean velocity (vi os ) is 
nowhere greater than that of the complete rotating model 
in which all stars rotate in the same sense. 



4 MODELLING M32 

Because of its high central density and steep central gra- 
dients in rotation and velocity dispersion, the nearby com- 
pact elliptical M32 is believed to lodge a massive black hole 
(Tonry 1987; Dressier & Richstone 1988; Richstone, Bower 
& Dressier 1990). However, this presumption was based on 
spherical mass models, despite M32 being as flat as E3. Re- 
cently, high spatial resolution data have been obtained both 
for the surface brightness by HST observations (Lauer et al. 
1992) and for the kinematics (VPs) by ground based spec- 
troscopy (van der Marel et al. 1994a, hereafter vdM94a). 
These data make it worthwhile to investigate the dynamics 
of M32 with self-consistent distribution-function models, as 
described above. Very recently, van der Marel et al. (1994b) 
have used a simple stellar mass model fitting the HST data 
to solve the Jeans equations up to third order assuming a 
2I-DF, while Qian et al. (1994) computed the actual DF for 
that model using the Hunter & Qian (1993) technique. I 
discuss the relation to this work in 4.6. 

4.1 The luminosity density 

Figure 2 shows the surface brightness and ellipticity profiles 
from which an axisymmetric luminosity density was derived. 
Fig. 2 combines the HST observations of Lauer et al. (1992) 
for the inner 4" (using V-R = 0.4 mag), the data of Peletier 
(1993) out to 32", and Kent's (1987) data further out. Inside 

1 This also holds for axisymmetric three-integral models provided 
the third integral is even in v^. 

2 An exception is the case in which the dominating mass compo- 
nent is natter than the stellar distribution, which reduces the need 
of near-circular orbits to flatten the stellar system, and hence in 
projection results in less motion on the major axis. 



Table 1. Deprojection of the surface brightness of M32 

* E 7 P Wrl XI [mag] 

90° E2.7 1.53 3.3 4 0.0057 
55° E4.5 1.53 3.3 5 0.0078 

The first two columns give inclination and intrinsic ellipticity; 
7 and f3 are the parameters of the first guess for the density 
(equation 6). -/Vrl is the number of RL iterations, while xi gives 
the RMS relative deviation in surface brightness. 

Table 2. The models: mass model and even part of the DF 



model 


i 


Th 


M. 


JVrl 


X P 


A 


90° 


2.4 




18 


0.0067 


B 


55° 


2.6 




18 


0.0065 


C 


55° 


2.6 


1.63X10 6 


21 


0.0062 


D 


55° 


2.6 


1.95X10 6 


21 


0.0063 



The first column labels the models; i is the inclination angle, T^ 
and M 9 are the stellar mass-to-light ratio in R and the mass of 
the central black hole, respectively, both in solar units. -/Vrl is 
the number of RL iterations performed, while \p gives the RMS 
relative deviation in the density. 

0'.'l5 a power law cusp Jocr -0 ' 53 with axis ratio 0.73 is as- 
sumed, which is consistent with the HST data (Lauer et al. 
1992). The deprojection has been carried out as explained in 

3.1 with the isophotes assumed to be exactly elliptical. Two 
angles of inclination have been considered: i = 90° (edge- 
on) and i = 55°. The details of the deprojections are given 
in Table 1. By more iterations the deviations from the in- 
put luminosity could be further reduced, however, this does 
not seem wise since the measured luminosity is hardly more 
accurate than 0.01 magnitudes. 

4.2 Constructing f e (E, L z ) 

In the case of i = 55° f e was constructed for three values of 
M, , the mass of the central black hole: 0, 1.63, and 1.95 mil- 
lion solar masses, while for i = 90° the galaxy was modelled 
with M % = only (see Table 2). All models have constant 
stellar mass-to-light ratio. The differences between p gener- 
ated by the DF and that obtained from the deprojection are 
of the order of (or smaller than) the uncertainty of the den- 
sity due to incomplete deprojection and erroneous surface 
brightness (see above). 

For models B and C (Table 2) Figure 3 shows f e as 
a function of Re for various fixed values of L 2 /L 2 ilc . The 
asymptotic behaviour with energy at small and large radii 
can easily be understood. In the envelope a power law for 
the density of pocr -3 8 is adopted, resulting in a scaling of 
/ oc R~j^ 3 in the keplerian regime. In the centre the density 
scales as r -1 ' 53 , which in the case of a black hole (keplerian 
again) gives / oc i?^ 03 , while for a self-consistent cusp one 
finds focR~ 2 235 (Dehnen 1993). 

The independence of the DFs can be understood by 
locally comparing the flattening of density and potential: 
the flatter the density compared to the potential and the 
steeper the density profile, the more high-angular momen- 
tum orbits must be occupied (Dehnen & Gerhard 1994; Qian 
et al. 1994). In the envelope of both models of Fig. 3 the po- 
tential is dominated by its monopole and almost round while 
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Figure 3. Even part of the 2I-DFs for models B and C plotted 
vs. R E , defined by y(R E ,0) = E, for L 2 z /L 2 chc = 0, 0.25, 0.5, 
0.75, 1 (from bottom to top). 

the density is moderately flattened (axis ratio PS 0.81) and 
steeply falling off. In the centre of model B (no black hole) 
the potential is flattened and the density is rather flat (axis 
ratio PS 0.55) but has a shallower gradient than in the en- 
velope. All together lead to a weaker independence in the 
inner parts of model B than in its outer parts. The opposite 
is true for model C with a central black hole, since the lat- 
ter makes 1 F spherical in the inner parsec and makes many 
high-i z orbits necessary. 



The DFs resemble each other in their small-scale be- 
haviour, which depends on the small-scale behaviour of 
p(^,R). There are two clear signatures. First, the wiggle 
at Re ~ 1 pc corresponds to the turn-over of density and 
surface brightness to a power law inside 0'.'l5. This gives 
a drastic change in dp/d^fr, what is closely related to the 
DF's run with energy, see Hunter & Qian's (1993) formula 
for f(E,L z ). A central power law in the density may have 
evolved during the adiabatic grow of a black hole (Young 
1980; Quinlan, Hernquist & Sigurdsson 1994). Therefore, the 
coincidence of the corresponding pattern in the DF with the 
energy, at which the DF changes its slope due to the pres- 
ence of a black hole, may be not accidental (note that M, 
is solely fitted to the kinematics, see below). However, in 
case of a adiabatically grown black hole, the DF should turn 
smoothly from one regime into the other, implying that the 
radial turn-over in the mass models used here is somewhat 
too drastic. Second, the wiggle at Re ~ 20 pc, which only 
occurs at L z ss L clvc , is related to the sudden increase in 
ellipticity at ~ 4" (Fig. 2), and hence may be somewhat 
artificial. 



4.3 The observables of f e (E, L z ) 

VdM94a have measured the VPs for several slit positions and 
fitted a Gauss-Hermite series to each of them (van der Marel 
& Franx 1993), from which reliable estimates for the RMS- 
velocities and the even parts of the velocity profiles, VP e , 
can be achieved. Fig. 4 compares the RMS-velocities with 
those of the constructed models, while Fig. 5 compares the 
VP e at some selected position on the sky. They have been 
computed using a Monte-Carlo technique, which takes into 
account effects of binning and seeing (Appendix C). The 
(constant) mass-to-light ratios T_r of the different models 
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Figure 4. RMS- velocities of M32 on four slit positions as observed 
by vdM94a and of the 2I-models (see Table 2). 



have been fixed to match the data on the minor axis. The 
following points regard further attention: 

First, the edge-on model (A) cannot match the kinemat- 
ics of M32: its VP e at 5" on the major axis is insufficiently 
fiat-topped (Fig. 5), and on the major axis it shows less 
RMS-motion than M32 (Fig. 4). The latter is untypical for 
elliptical galaxies, which usually have a somewhat greater 



ratio 



\l/2 



l*Tos)min/( ? 'ios}maj than an 2I-model seen edge-on (van 
der Marel 1991). The models with i = 55° show both cor- 
rectly the ratio of minor to major axis motion and the shape 
of VP e outside the centre. 

Second, the models without 
to reproduce the high velocities 
even have centrally declining 
consistent stellar cusps (Binney 1980; Dehnen 1993). On the 
other hand, the models with central point masses can ac- 
count for the measured high central (ffos) 1 ^ 2 (in particular, 
the central VP is in excellent agreement with model D - 
see Fig. 5). However, they have slightly more RMS-motion 
at projected radii of 1-2". This is clearest on the minor axis 
but not very significant if one compares the VP e in Fig. 5. 
A 2I-model with a less massive black hole would have the 
correct (ffos) 1 ^ 2 at 1-2", but would be unable to account for 
the high observed velocities in the very centre. Hence, two- 
integral models cannot reproduce this effect, which therefore 
is a hint that the DF depends on the third integral. 

Third, at about 3-5" on the minor and intermediate 



a black hole are unable 
inside l" (Fig. 4); they 
as is typical for self- 



vl/2 




Figure 5. Comparing the even parts of the VPs (normalized to 
unit area) in the centre, and at l" and 5" on the major, in- 
termediate and minor axes. The VPs of M32 are represented by 
points, which have been drawn from the the Gauss-Hermite fit 
by vdM94a assuming that the given errors are uncorrelated, nor- 
mally distributed, and are the only errors. In the upper three 
panels most of these points are overlayed by the VPs of model C 
or D. 



axes the models show less RMS-velocity than observed. How- 
ever, in the VP e (Fig. 5) this appears to be only a marginal 
effect, and does not occur on the major axis, where the errors 
are smallest. 



4.4 The observables of complete distribution 
functions 

Instead of trying to find the odd part for the DF that results 
in the best agreement of the predicted and observed VPs, I 
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Figure 6. Parameters of Gauss-Hermite fits to the VPs of M32 measured by vdM94a. From top to bottom: mean and dispersion of the 
best-fitting Gaussian to the VP, and the Gauss-Hermite coefficients up to order six. The curves show the corresponding values of the 
models C (solid) and D (short dashed), which lodge black holes of masses 1.63 X 10 6 Mq and 1.95 X 10 6 Mq, respectively, and of model 
B (long dashed) without black hole. The atmospheric seeing as well as the spatial binning of the data have been simulated. 



make the simple ansatz 



f(E,L z ): 



= /e(£,0) 



L z sC 



2f e {E,L 2 z )-e aL '' L °"°f e {E,0) L z > 0, 



: «i + (a 2 — «i) 



R% + Rl 



(9) 



(10) 



where a\, ct2 and Ro are free parameters. This ansatz im- 
plies, that the DF falls off exponentially with L z /L clvc for 
retrograde orbits with a scale a -1 that depends on energy 
in a very simple form: a = a\ for Re <C -Ro and a = ct2 for 
Re 3> -Ro with a smooth transition. After a few experiments, 
I found the models with a\ = 6, a\ = 1.7, and Ro = 8 pc 
to give kinematics in reasonable agreement with the data of 
vdM94a. Figure 6 compares the Gauss-Hermite-fit parame- 
ters (van der Marel & Franx 1993) vm, am and hs, . . . ,h% 
measured for M32 by vdM94a with those of the models B, C 
and D, while Figure 7 shows the VPs along the major axis of 



the models C and D (those with a central black hole) and of 
M32, as reconstructed from the Gauss-Hermite-fit parame- 
ters and their formal errors. 

Even for the simple form of f given by equation (9), 
the agreement of the models containing a black hole with 
the VPs of M32 is remarkable, much better than one would 
expect based on the large variety of three-integral distribu- 
tion functions and their observables, that are possible in an 
oblate stellar system (Dehnen & Gerhard 1993). The follo- 
wing points are worth noting. 

First, beside leading to the high central rotation and 
dispersion, the black hole also influences the VP-shapes, 
which is obvious especially from the /t 8 -profiles on the ma- 
jor axis: the dips in A3 at l" and in A4 in the very centre 
are not present in the model without central point mass. 
The negative hi for the central VP is not what one would 
naively expect, since the fast motions near the black hole 
should create broad high-velocity wings resulting in positive 
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Figure 7. Comparing the VPs (normalized to unit area) along the major axis of the models with central black hole with those of M32 
as observed by vdM94a. The latter are represented by points as in Fig. 5. 



hi. However, seeing smears in the VPs at ~± 0'.'7, which 
show rotation in opposite directions and thus broaden the 
low- velocity region of the central VP leading to negative hi. 
Higher spatial resolution would yield positive hi for the cen- 
tral VP, see subsection 4.5 below. 

Second, the predicted vm, crm and h t differ slightly from 
the observed values, which is clearest on the major axis, 
where the errors are smallest. However, this seems to be 
hardly significant from the VPs in Fig. 7, where only the 
formal errors are considered and no systematic errors (e.g. 
template mismatching). Additionally, the VP-shape parame- 
ters h t are quite sensitive to f , as experiments with various 
a (equation 10) have shown, indicating that there might well 
be a f which fits the VPs even better. It is interesting to 
note that the better fitting models have relatively less retro- 
grade stars in the centre than outside about 7-10 pcw2-3", 



i.e. rotate relatively more strongly in the region where the 
black hole reigns. 

Third, there are few minor discrepancies, which may in- 
dicate that f depends on ^3, in particular the rotation at 
slit position 4" parallel to the major axis. The non-zero val- 
ues of A3 and hs on the minor axis, as well as the different 
rotation velocities at the same central distances on the two 
intermediate axes are even not consistent with any axisym- 
metric equilibrium model. However, it is not clear that these 
discrepancies are statistically significant. 

The main conclusion is that a 2I-model with i PS 55° and 
M. = 1.6-2 x 10 6 Mq is mainly consistent with the kinematic 
data. This makes it unlikely, that any 3I-DF which fits the 
data equally well depends strongly on the third integral. 
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Figure 8. Predictions of the models with central black hole for 
observations with the r /09 X r /09 aperture of the Faint Object 
Spectrograph of the HST: the Gauss-Hermite fit parameters of 
the VPs along the innermost arcsecond of the major axis. 



4.5 Predictions for observations with the HST 

As has already been pointed out by Qian et al. (1994) HST- 
spectroscopy of the centre of M32 would be worthwhile, es- 
pecially since the integration times necessary to gain high 
signal-to-noise ratios are short because of the high surface 
brightness. The models can be compared with such obser- 
vations to further constrain the mass of a massive central 
object in M32 (note that the O'.'l resolution of the HST cor- 
responds to PS 0.3 pc at the distance of M32). 

For models C and D (with black holes) I have com- 
puted the kinematics on the major axis inside l", which 
would be measured using a square aperture of 0'.'09 x 0'.'09, 
the smallest available aperture on the HST Faint Object 
Spectrogaph. Figure 8 shows the corresponding values of 
the Gauss-Hermite-fit parameters along the major axis. At 
such a high resolution the central black hole would leave 
clear imprints on the kinematics: a high central velocity dis- 
persion and strong gradient in projected rotation, as well as 
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Figure 9. The major axis' VPs (normalized to unit area) at 0" 
and 0'/2 that the models with central black hole predict for obser- 
vations with the 0'/09 X 0'/09 aperture of the Faint Object Spec- 
trograph aboard the HST. 

a central VP showing strong deviations from gaussian form 
- see also Fig. 9, where the predicted VPs at 0" and 0'.'2 
are plotted. The high velocities in the neighbourhood of the 
black hole lead to broad wings of the central VP measured 
by values of hi > 0.05 (in contrast to the earth-based obser- 
vations of vdM94a, for which hi < in the centre, see the 
discussion above). A measurement of central VPs like those 
predicted here would further constrain the possible stellar 
dynamical configurations and provide more evidence for the 
presence of a black hole in M32. Another clue for a central 
dark mass would be the detection of stellar velocities which 
are well above the central escape velocity v* sc of the po- 
tential due to the stars alone. Such velocities would imply 
a deeper central potential and hence a dark component at 
least as compact as the stellar distribution - they cannot 
be explained by an extended dark halo, since then the as- 
sociated low-binding-energy stars should also be visible at 
much larger radii. However, for M32 such a detection seems 
to be hardly possible, since the predicted fraction of stars 
with projected velocities exceeding v* sc PS 387kms -1 is very 
small, see Fig. 9. 

4.6 Comparison with other recent work 

Van der Marel et al. (1994b) used the mass model 

-1.435 /-, , r / n ;; cc n2\-0.423 /-,-,-, 

poem (1 + [ra/0. 55J ) (11) 

with m 2 = R 2 +(2/0. 73) 2 (assuming edge-on projection) to 
solve the Jeans equations up to third order, where a cen- 
tral black hole with mass 1.8 X 10 6 Mq was added to the 
potential. Later on, Qian et al. (1994) have used the Hunter 
& Qian (1993) method to evaluate the 2I-DF for this mo- 
del, and subsequently computed the VPs (which were not 
available from the pure Jeans equation model). This model 
differs from the successful models C and D above mainly in 
the four following points. It (i) assumes edge-on inclination, 
(ii) has a somewhat shallower cusp, (iii) has constant ellip- 
ticity, while M32 becomes rounder than E2.7 already at 10", 
and (iv) its surface brightness declines as r -1 ' 28 outside a 
few arcseconds, while that of M32 falls off steeper outside 
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about 20". Nevertheless, their results are very similar, even 
in some details, to those found in this work. However, there 
are two noticeable differences. 

First, as for the edge-on model A above, they found 
too much motion on the minor axis and VP e on the major 
axis being less flat-topped than those observed. The authors 
interpret the latter as a suggestion that f e depends on ^3, 
while in the models presented here these discrepancies were 
removed by taking i PS 55°. However Qian et al. claim this 
discrepancy in the VP e is hardly significant, in apparent con- 
flict with Fig. 5. 

The A 8 -profiles of the major axis differ: A3 of the Qian 
et al. model does not show a dip at f ", as do the models pre- 
sented here, and A4 is smaller at all radii. Both differences 
may well be caused by the different forms for f that have 
been used. However, in the centre there are very few retro- 
grade stars in either model, and hence the freedom in f is 
small. It seems possible that the discrepancy in A3 is due to 
the shallower cusp of their model, since the steeper the den- 
sity profile the more high-i z orbits a 2I-DF has to populate, 
and hence the more asymmetric the intrinsic and projected 
velocity distributions will be, leading to more negative A3. 



5 SUMMARY AND CONCLUSIONS 

The Richardson (1972) -Lucy (1974) algorithm can be used 
to quickly and accurately compute the two-integral distribu- 
tion function f(E, L z ) for an axisymmetric stellar system; 
there are no restrictions regarding the functional form of the 
density p(R, z) or otherwise, as long as density and poten- 
tial can be represented on a grid, as is the DF itself. The 
computed 2I-DF is by construction positive and hence phy- 
sical. Nevertheless, even unphysical DFs can be obtained as 
the difference of two positive ones. For better convergence of 
the iterative RL algorithm, it is very helpful, to transform 
the integral equation to be solved in an optimum form, in 
which the known projected quantity (e.g. the density) is as 
homogeneous as possible, and the kernel (e.g. the density 
of the orbits with same E and L z ) is as local as possible 
(Appendix A). 

In Section 3 it is shown, how this technique can be used 
to construct 21- DFs for individual galaxies with the observed 
surface brightness distribution as the only input. Because 
of the uniqueness of f e and the simplicity of its evaluation 
(compared to 31- DFs), the 2I-model can be used as a tem- 
plate model for the comparison of the kinematic data of 
the galaxy under investigation. In particular, it should be 
possible to tell whether the velocity distribution is less or 
more tangentially or radially anisotropic compared to that 
of the 2I-model, and if so, how strong the effects are, i.e. 
how strongly the DF depends on ^3 . 

In Section 4 the technique has been applied to the 
nearby elliptical galaxy M32. Starting from a composed sur- 
face brightness, which comprises the data of Kent (1987) 
and Peletier (1993) for r > 4", the HST data of Lauer et al. 
(1992) for O'.'f 5 < r < 4" and a cusp I oc r"° 53 for r < O'.'f 5 
as suggested by Lauer et al., the axisymmetric density was 
obtained for the inclination angles 55° and 90° (edge-on). 
For each of the resulting stellar distributions, I have first 
computed the 2I-DF without assuming a central black hole 
and have evaluated the projected kinematics (VPs), as if 



obtained under the observing conditions of vdM94a. The 
edge-on model was unable to fit the observed kinematics: 
it shows too little motion on the major axis in comparison 
to the minor axis, and significantly different shapes for the 
VPs on the major axis. The 55° inclined model (intrinsic 
E4.5) gave both correctly, but was unable to reproduce the 
high observed velocities inside 2". To account for the latter, 
two additional models have been constructed with central 
point masses of 1.63 x 1O 6 M and 1.95 x 1O 6 M . For the 
odd parts of the DFs a simple ansatz has been made. These 
models have kinematic observables, including the shapes of 
the VPs, in excellent agreement with those of M32 given by 
vdM94a. This means that the available kinematic data are 
consistent with a 2I-model with central black hole (or other 
compact dark mass) in the mass range from 1.6 to 2 million 
solar masses. It is not clear, whether three-integral models 
exist that equally well fit the observed kinematics, but con- 
tain a less or more massive black hole. Nevertheless, it seems 
unlikely, that such models would manage without any cen- 
tral dark point mass, for the following reason. To explain the 
high central dispersion and rotation without black hole, a 
highly radially anisotropic DF is required with the high cen- 
tral velocities occurring at the pericentres of stars on highly 
eccentric orbits with low binding energy. On the other hand, 
the kinematics of M32 outside of the innermost arcseconds 
indicate tangential anisotropy there (since a E4.5 2I-model 
fits the data), which is certainly not consistent with many 
stars being on highly eccentric orbits. From this argument I 
infer the presence of a dark mass of 1-3 x 10 6 Mq in the centre 
of M32. Kinematic data at the resolution of the HST would 
give further constrains on the possible DFs of M32. For ob- 
servations with the 0'.'09 x 0'.'09 aperture of the HST Faint 
Object Spectrograph the 2I-models with black hole predict a 
clearly non-gaussian central VP with broad wings, true and 
gauss-fit velocity dispersion of 151 kms -1 and 121 kms -1 , 
respectively, for a black hole mass of 1.63 x 10 6 Mq, while 
a black hole of 1.95 x 1O 6 M would give 171 kms -1 and 
131 kms" 1 . 

The possibility to model not only the mean velocity and 
dispersion of the line-of-sight velocity distributions of gala- 
xies, but also their shapes, the velocity profiles, is a great 
challenge. The additional information in the VPs constrains 
the possible dynamical configurations of the observed ga- 
laxies - as this paper demonstrates for M32 - and hence 
will lead to a better understanding of their internal struc- 
ture. The main problem is the construction of dynamical 
models sophisticated enough to exploit the information in 
the VPs. The simplest such models are axisymmetric with 
/ = f(E,L z ). Only for very few galaxies have distribution 
functions been evaluated so far, and further work in this 
direction is certainly needed. 

I anticipate making my FORTRAN codes for the evalua- 
tion of f e (E, L 2 Z ) and of its VPs available to the community 
within a few month after publication of this paper. 
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APPENDIX A: THE MODIFIED 
RICHARDSON-LUCY ALGORITHM 

Al The original algorithm 

Richardson (1972) and Lucy (1974) described an algorithm 
(hereafter RL algorithm) to solve the projection equation 



for the unknown (intrinsic) ip(£.) given the known projected 
(observed) function <p(x) and the kernel K(x\£), which all 
have to be non-negative and normalizable. Here the co- 
ordinates x and { may well be multi-dimensional. The 
scheme consists of successive applications of the two iter- 
ative steps: 



<Pn(x) 



Mt)K(x\t)dt 



f[<p{x)/<Pn{x)]K(x\£)dx 
J K(x\£)dx 



(A2) 



(A3) 



where the subscript n indexes the iteration step. In equation 
(A3) the integral in the observed domain includes all points 
x to which ip(£.) contributes. This algorithm converges in 
the sense that the logarithmic likelihood of the projected 
function, 



tp(x) log[tp n (x)]dx, 



(A4) 



is maximized under the constraint J ip(x)dx = J ip n (x)dx 
(Shepp & Vardi 1982). 



A2 Modifying the RL algorithm 

The RL algorithm may easily be modified by manipulating 
the projection equation (Al) such that its solution remains 
unchanged. The simplest way to do so is by multiplying 
it on both sides by some positive function a(x), which is 
only restricted by the constraint, that the product a(x)ip(x) 
must remain normalizable 3 . This manipulation affects only 
the correction step (equation A3), which becomes 

I \ip(x) /tpn(x)] a(x) K(x\£)dx 

f n+1 (o = Mo n 1 r \L\L A • ( A5 ) 

J a(x) A (x\£)dx 
The such modified algorithm converges in the sense that 



H n (a) 



a(x) ip(x) ln[ip n (x)]dx 



(A6) 



ip{x) 



W)K(x\t)dt 



(Al) 



is maximized under the constraint J ip(x)dx = J ip n (x)dx. 

Such a modification is sensible if the given projected 
function <p(x) is significantly inhomogeneous, since in this 
case H n (equation A4) is dominated by points where <p(x) 
is greatest, with the result that in regions where <p(x) is 
small the original RL algorithm tends to converge rather 
slowly. If, for instance, <p(x) is nowhere zero, a suitable choice 
is a(x) = ip(x)~ 1 , which results in a scheme maximizing 
j dx ln[tp n (x)/tp(x)]. 

Another problem often encountered with the RL algo- 
rithm is the fitting of noise, which typically occurs after 
the first few iterations. Numerically, the unknown ip to be 
found is always represented on some grid, and hence only 
given with finite resolution. Therefore, one may smooth ip 



3 Another method would be to differentiate the projection equa- 
tion with respect to some components of x (which in the simplest 
case ip = Tpd£ already gives the solution) or any combination of 
this with multiplication. However, this method is useful for the RL 
algorithm only if the left hand side of the new projection equation 
is everywhere non-negative. In practice the accurate evaluation of 
the involved derivatives is likely to be problematic. 
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after each correction and before the projection step with a 
smoothing scale similar to the resolution element of the grid. 
This will guarantee a smooth result and avoid the fitting of 
small-scale noise. 



A3 Accelerating the RL algorithm 

For the purpose of speeding up the convergence, at each it- 
eration the new guess for the intrinsic function ip n +i(£,) is 
written as ipn(£.) + 8,p(^), where 8,p(^) is the correction of 
ipn(^), which projects into 8 <p (x) = J 8^ ({) K(x\^) d{. For 
accelerating the convergence of the algorithm, 8,p(^) is mul- 
tiplied by a factor v , which after projection is chosen to be 
optimal and to obey 1 ^ v < /i/ mlI , where J'max is given by 
the condition that ip n +i(£,) must not become zero and / ^ 1 
is a fudge factor. 

Usually, the optimal v is determined by maximizing the 
resulting _ff n +i (Kaufman 1987; Lucy 1992). This can be 
done by numerically solving &H n -yi / &v = 0, which however, 
again involves an iterative procedure. Another possibility is, 
to choose the optimum v by minimizing the relative RMS- 
error x between tp and <p n +i, which yields 



I'opt 



(f ~ Pr, 



E 



(A7) 



where the sums go over all data points. 

Clearly, the function ipn~ which maximizes the likelihood 
H will be distinct from the function ip x which minimizes the 
relative RMS-error, and therefore the choice of v as in equa- 
tion (A7) may retard or even stop the convergence. However, 
for the first iteration steps ip n will differ clearly from both 
functions ipn~ and ip x when compared to their difference, and 
hence maximizing H or minimizing \ wm be quite similar. I 
found this choice of v op t to work quite well, especially if the 
RL algorithm is suitably modified, which gives ipn~ ~ ipx- 



APPENDIX B: ALTERNATIVE RL SCHEMES 
TO RECOVER f e (E,L 2 z ) 

As mentioned in Section 2, one can apply the RL algorithm 
on the basic integral equation to be solved after taking a 
derivative on both of its sides, which reduces the dimension- 
ality of the involved integrals by one, and hence should result 
in numerically faster RL schemes to recover f e (E, L 2 ) than 
that given in Section 2. However, the accurate knowledge of 
the derivatives involved is required. 

Bl Scheme 1: taking the derivative 9/9\I/ 

Taking the partial derivative with respect to 'F on both sides 
of equation (1) results in the integral equation 



A(^,R) 



Qp(^,R) _ 9 3/2 



3* " "Jo VW^E 
The modified RL algorithm gives for the correction step 
r* 2 _c 



fe,n+l(E, L z ) 
fe,n(E,Ll) 

for L z ^0, and 



-E [ A^.R) j^^—^ 



(B.2) 



/e,„+l(£,0) /j 



*(0,0) d ,j 



A(y,oy- p /A n (y,o) 



fe,n(E,0) 



r*(0,0) d<t 



A(V,0)-p 



(B.3) 



for L z = 0. $i , ^2 are given in Section 2, while the exponent 
p can be chosen to optimize the convergence. 



B2 Scheme 2: taking the derivative 3/3i? 

Multiplying equation (1) with R and subsequently taking 
the partial derivative with respect to R on both sides yields 



B(f,R) 



d(Rp[$,R]) 
dR 



:2 5/ V 



c dE V* - E f e (E,2R 2 [^ - E]). (B.4) 
Jo 

The correction step of the modified RL algorithm reads 

r* 2 ^ [flft.fl) 1 -'" 

J*! u * I B n (t,R) 



R=L z l^/2(^-E) 



fe,n+l (E, L z ) 
fe,n(E,L 2 z ) 

for L z ^ and is given by equation (3) for L z = 0. 



f* 2 d* —r^= 

J *l 1 ' ,/2(¥- 



(B.5) 



APPENDIX C: 
PROFILES 



COMPUTING VELOCITY 



The velocity profile (VP) is given by a three-dimensional 
integral over the distribution function, namely over the two 
velocity components in the plane of sky and along the line 
of sight. However, if one wants to compare with observed 
VPs, effects of binning and seeing have to be included, each 
of them resulting in another two-dimensional integration. 
The best way of handling such high dimensional integrals is 
a Monte Carlo integration, a straightforward application of 
which gives a highly ineffective procedure for the evaluation 
of a VP, since most of the randomly chosen phase space 
points lie outside the region which dominates the VP, i.e. at 
small or even negative binding energies. 

To get a more efficient procedure one needs to preferen- 
tially pick out phase space points resulting in large values of 
the DF. In Monte Carlo integration technique this is called 
"reduction of variance" , it is equivalent to a substitution of 
a variable of integration, i.e. one has to multiply with the 
corresponding Jacobian. After some experiments I ended up 
with the iteration of the following steps: 

(1) Set F = 1, the factor containing the Jacobians due to 
the reductions of variance and the DF. 

(2) Choose a position on the sky (x,y) at random out of 
the binning area, and add a two dimensional vector 
drawn randomly from the point spread function. Both 
times one may well prefer small radii to account for the 
gradient in surface brightness, in which case F has to 
be multiplied by the corresponding reduction factors. 

(3) Choose a position s on the line of sight at random out 
of ^ s/(s + a) < 1; a = \/ x 2 + y 2 + a 2 , where a is 
a small number, and s = corresponds to the point, 
where a pure spheroidal density distribution with simi- 
lar ellipticity would be maximal along the line of sight; 
compute the galaxy's intrinsic co-ordinates (R,z). To 
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correct for the preference of small s, the factor F has 
to be multiplied by (s + a) 2 /a. 

(4) Choose energy randomly out of < E 2 ^ *(_R, z f , 
calculate v = ^2(* - E), and multiply F by 
■$(R,z) 2 /(2E). 

(5) Choose v,f, randomly out of —v ^ v$ ^ v, calculate v m = 
\Jv 2 — v 2 ^, evaluate f(E,Rv<p), and multiply F with 
vf{E,Rvt) 

(6) Loop over ^ ip < 27T, for each ip set «n = v m cosip, 
v z = v m smip, compute vi os , and add F to the corre- 
sponding pixel of the VP-histogram. 

Points (4) to (6) are based on the identity 

f /"* 2 r\F 2 /•V 2 (*-- B ) r 2n 

/ d3 "= / / , d ^ / 

A 2 ^* Jo ^ J-y/2(y-E) Jo 

The most time consuming piece in this whole procedure is 
the evaluation of the DF, since it involves a two-dimensional 
interpolation. The trick in the above algorithm is the loop in 
step (6), which exploits the fact, that a 2I-DF is independent 
on the angle ip in velocity space. This allows the number of 
evaluations of f(E, L z ) to be much smaller than the number 
of points actually added to the histogram. 
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